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We study gravitational perturbations of Schwarzschild spacetime by solving a hyperboloidal initial 
value problem for the Bardeen-Press equation. Compactification along hyperboloidal surfaces in a 
scri-fixing gauge allows us to have access to the gravitational waveform at null infinity in a general 
, setup. We argue that this hyperboloidal approach leads to a more accurate and efficient calculation of 

the radiation signal than the common approach where a timelike outer boundary is introduced. The 
method can be generalized to study perturbations of Kerr spacetime using the Teukolsky equation. 
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I. INTRODUCTION 



• In the last few years, numerical relativity has made tremendous progress in treating the inspiral of compact objects, 
q-i starting with the breakthrough simulations of P, 0, Q ■ Rather accurate simulations of the last ~ 10 orbits of roughly 
equal mass binary black holes have become standard in the field (see e.g. [1, [H, H, 0, H, Q for recent examples). 
However, these simulations are still restricted to finite outer boundaries and gravitational radiation is extracted at 
a number of finite shells. Maybe not surpris ingly, it has been found that the wave extraction radius is a significant 
limiting factor for the precision of results [lfj, |ll|, QJ, [l3| . The practical implementation of algorithms that include null 
infinity in the numerical grid is thus highly desirable. It also provides for a fruitful problem to foster the interaction 
' of numerical and mathematical relativity. 
£SJ , An important stepping stone to develop such algorithms is the case of linearized perturbations on a fixed background. 
0^ ' The numerical study of linearized perturbations is convenient as one can evolve the gravitational radiation very 
accurately for a comparatively low computational cost llil. [TH [Hj|. For evolutions along null surfaces this has 
. provided important guidance in the past, see e.g. (TtI [l8l.Tl3| . Note, however, that linear analysis may not be valid 
t-H ' in certain special weak field situations. Examples where the linear analysis breaks down have been found in studies 
of test fields [2(| In the pure grav itational case, so far it seems that the perturbative analysis fits well with the 
O ■ nonlinear evolution [H, S3, Hi, [Mill • 

When the background spacetime describes the particularly interesting case of a stationary black hole, the perturba- 
tion equations take the form of linear wave equations with a potential term (in our case the Bardeen-Press equation 
[11). The solutions to such equations show three stages: an initial transient (the "direct signal"), quasi-normal mode 



ringing and power-law decay. The standard method to numerically analyze these stages in the time domain has been 
to solve initial boundary value problems using Cauchy-type foliations. This method has certain disadvantages. When 
one is primarily interested in radiative properties of the solution far away from the source, the use of Cauchy-type 
foliations is a waste of computational resources. Furthermore, there arc difficulties related to the choice of boundary 
conditions. It has been shown that a typical choice of boundary data for numerical relativity destroys the polynomial 
tail behavior of solutions to wave equations on a Schwarzschild spacetime [13, HI . Note that even if the polynomial 
tail can be calculated by using a better choice of boundary data, the decay rate will be the one near time-like infinity. 
It has been pointed out, however, that the relevant decay rate for an astrophysical observer is the one along null 
infinity [H, M, [p • 

Using foliations that approach null infinity would cure these problems. Numerical studies of gravitational per- 
turbations at null infinity have been performed by [tR [l8| using a characteristic approach for the Bardeen-Press 
equation, and have been compared with fully nonlinear studies in [jjjj . The authors have been able to efficiently 
calculate the quasi-normal mode ringing and the power law decay. There are, however, two main difficulties in the 
characteristic approach. An important aspect of solving the Bardeen-Press equation instead of the Regge-Wheeler- 
Zerilli equations is that the favorable separabil ity pro perties of the Teukolsky formalism extends to Kerr spacetime 
[32| . Unfortunately, generalizing the studies of [l7l. Il8l| to Kerr spacetime is complicated because the construction of 
null foliations smoothly covering Kerr spacetime is not straightforward and from a numerical point of view not very 
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convenient [33l. |34|. A further complication is related to the causal structure of a black hole spacetime. Due to the 
global relationship between null surfaces on which boundary information is known to pose a characteristic initial value 
problem, the numerical solution for a perturbed black hole needs to be calculated in two steps: advanced [HI and 
retarded [l7| ■ A related difficulty shows itself also in the choice of suitable coordinates which allow the computation 
of the gravitational waveform at null infinity and provide good resolution near the horizon. 

An alternative, and very flexible method to study radiation is to use spacelikc surfaces that approach null infinity. 
Following [35l ] we call such surfaces hyperboloidal and the related initial value problem a hyperboloidal initial value 
problem. Hyperboloidal surfaces are as flexible as Cauchy-typc surfaces because the property that characterizes them 
is not local but global. Further, they enable, like characteristic surfaces, a clean treatment of gravitational radiation 
because they approach null infinity. In this sense, they are hoped to combine the "best of both worlds" , and we expect 
increasing interest in this approach. 

In the present paper we want to learn from the linearized case as a preparation for the nonlinear hyperboloidal initial 
value problem. To this end we construct solutions of the Bardeen-Press equation by numerically solving a hyperboloidal 
initial value problem for the perturbations of Schwarzschild spacetime and we calculate the gravitational wave signal 
at null infinity. We find that compactification along hyperboloidal surfaces leads to a very efficient code. Specifically, 
we show that the calculation of the quasi-normal mode frequencies in the time domain can be performed accurately 
with a small number of grid points using standard numerical techniques. The last stage of the solution, namely the 
power-law decay, can be numerically calculated if the resolution is sufficiently high. A major advantage of our method 
for this part of the signal is not so much the efficiency but rather the fact that we can measure the decay rate at null 
infinity by solving a Cauchy problem in the PDE sense. 

In the presentation and, to some extent, in the choice of the methods applied in the paper, we tried to avoid 
a heavy use of the conformal language. Our philosophy of constructing an algorithm for the hyperboloidal initial 
value problem follows the suggestions in [36L \3l\ . Instead of constructing a gauge- independent, regular, conformal 
Bardeen-Press equation, we make a suitable use of our gauge freedom before compactifying the coordinates and show 
that, in the end, regular equations are obtained. While both approaches are equivalent in their final aim, we hope 
that the way we present the method will be more accessible to numerical relativists who are not familiar with notions 
related to conformal geometry. 

The paper is organized as follows: In section II we present the Teukolsky formalism and derive the Bardeen-Press 
equation for a spherically symmetric but otherwise arbitrary foliation. In the third section we discuss our method 
to include null infinity on a numerical grid. We show that the hyperboloidal Bardeen-Press equation written in 
compactifying coordinates is regular at null infinity. After a short presentation of our numerical methods in section 
IV, we go over to the numerical studies of the last two stages for the gravitational perturbations where we demonstrate 
some of the basic advantages of the hyperboloidal approach in numerical work. In the discussion section we summarize 
our results and point out some ideas for future work. 



Teukolsky [32| derived a master perturbation equation for various fields in a Kerr background, using the spinor 
formalism introduced by Newman and Penrose [38| . One of the central ideas of this formulation consists in choosing 
a tetrad of null vectors, Z a ^ , and then defining the directional operators as the projections of the covariant derivative 
along each null vector. For a complete review on the subject and its properties, the reader is referred to [39L [40|. 

Usually the tetrad is defined with four null vectors: two real and along the light cone, the other two complex 
conjugate and in the perpendicular plane to the cone. There is much room for conventions regarding the orientation 
and scaling of the null vectors, and consequently for confusion. We make the following choices: we label the two real 
vectors along the light cone as = l^,Zi^ = fc M and choose them both future directed with ? M pointing outward 
and fc M pointing inward. The complex two null vectors will be denoted by m* 1 : Zi* 1 = m/', Z^ = m*^. 

After these choices, there is still the issue of the normalization of the products k^, and m M m* M . The rest of the 
products are zero by construction. The properties mentioned above can be expressed in the following equations: 



II. THE PERTURBATION EQUATION 



1. General definitions 



Z a ^ Zt>^ — rjat,, 



9^ = 



2rf 



a b 



z, 



(i) 



with r] a b a matrix of the form 
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with 7?i2 = ±1 depending on the signature. If we choose the signature of the spacetime to be (+,—,—,—) as in 
[1^, Ho|, then the choice consistent with Eqs. ([TJ) to derive the perturbation equation is 7y 12 = 1 [32l |. If the spacetime 
is described using the signature (— , +, +, +) we have to choose 7712 = —1 in order to remain consistent with Eqs. |T]). 
Some discussion on the change of the tetrad due to signature can be found in [4l| . 

We derived the linearized equation for the perturbed Weyl scalar "J 4 (note that the background value vanishes), 
leaving uns pec ified the choice of signature. Using the definitions for the directional operators and spinor coefficients 
as given in [39l . |40| , the final perturbation equation takes the form 

((A + r] 12 (4/i 5 + /j s * +37 s - 7s *)) (D - 7712 (p s -4e,)) - (5* + 7712 (3a s + [3 S * +4tt s - t s *)) (S + r) l2 (4ft - r s )) 

- 3 ? 712* 2 )*4 = (3) 

We have added a subindex s to the spinor coefficients in order to avoid confusion with other quantities that will be 
used later. In the rest of the work we will be using the (— , +, +, +) signature, and thus 7712 = — 1. 



A. Spherical symmetry 

We will work in spherically symmetric static spacetimes. Then the metric can be written as 

g = (-a 2 + 7 2 /3 2 ) dt 2 + 2 1 2 (3dtdr + 7 2 dr 2 + r 2 da 2 . (4) 

We will write the Bardeen-Press equation for \I/4 on such a background with respect to an orthonormal null tetrad 
(I, k, m, TTi), consistent with the conditions mentioned above. The tetrad can be given by 

*" = 2c?( 1, 7 _j8 ' 0,0 )' fcM " irf +l3 ) ,0 '°) ' "»" = -^-(0,0,l > icBcff) ) fh^ = {m ti )*. 

Remember that ^4 is defined by contracting the Weyl tensor with respect to the inward pointing null vector field of 
the null tetrad, k in our case. 

We compute the corresponding spinor coefficients and directional derivatives and substitute into the perturbation 
equation ([3]) . Using the definitions for the spinor coefficients given in (39j , and denoting derivatives with respect to r 
by a prime, we obtain the non zero spinors and Weyl scalars as 

Ps — ira 1 ' Ps — r j t s — , J S — 2q7 , 

o _ cote . _ a q, _ ( r2g 7fe M.)' /c) 

Before proceeding to write down the perturbation equation for ^4, we first note that the simple sourceless wave 
equation, 

g^V ll V^(x a ) = 0, (6) 
independently of the tensorial or spinorial character of <& in this background takes the form 

[Oh + O dv ] $ = 0, 



d 2 i a d 2 ( ( a\ 2 2 | ^ 2 a / r 2 /? 7 \ ' <9 a ( 2 ( a (3 2 j\\ d 



with 

rt a 2 ^ dt 2 ^ dt dr I ^ 7 ) ^ J dr 2 ~^ r 2 7 a ) dt r 2 7 ^ ^7 a ) ) dr J ' ^ 
d 2 1 d 2 n d 

° e * = W 2 + ^reW 2+cot6 de- (8) 

In the perturbation equation (|3|) we use the sign consistent with our choice of signature, that is 7712 = — 1. After 
substituting our tetrad and the obtained directional operators, spinor coefficients and the Weyl scalar ^2, we obtain 

[L tr + L ev ] #4 = 0, (9) 
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where the operators L tr and Lg^ are given by 



r 2 q2 r 2 q2 / a 2 _ p2^2\ q2 q g 

-J2 p»2 @ a 2 QiQ r y a 2 ) dr 2 dt~^~ dr 

cos 8 d ^ cos2 1 
sin 2 8 dip sin- 
Here we have defined 



Lev = Og v -4i — — Q — -2 . 



a 7 \ a J a d V r / 

B = (r ajpsHs) -4r — + 2r d - , 

«7 \ a / a \r J I 

C = — - — (3 {a"f p s ij, s r 2 Y + p s r (a ~/ fi s r)' — 10 fi s r (a 7 p s r)'^ — // s ( 5/> s + r/9 s ' +2r ( i^l^L ) ) 
r«7 V / \ y 07 y y 

2 

Q! 7^ 

There are several things worth noticing in this equation. Its principal part is the same as of the wave equation. The 
gravitational contribution is described by some extra coefficients in the terms containing no higher than first radial 
derivatives. The left hand side of the equation for the gravitational perturbation, as well as that for the wave equation, 
is separable into terms involving the (r, t) and the angular coordinates. The equation is somewhat lengthy, due to 
the fact that we are considering general metric coefficients, but nevertheless, it is quite tractable. Finally, we remark 
the fact that had we not considered the corresponding change in signature, we would have obtained an inconsistent 
expression. 

We will work with the function ^4 = r ^4, in order to factor out the expected decay due to the peeling behavior. 
The function ^4 satisfies the equation: 



L fr + Lf) 



*4 = 0, (13) 



where the operator L tr reads 



r 2 d 2 . „ „ r 2 d 2 . „ 4 d 2 ~ d ~ d ■ 



L *r = 2^2+ 2 ( 3 —7^7r + 2r Psf I s 7 ^+A- + B—-2r z C, (14) 

or dt or at or dr z dt or 



where now 



A = A-2r^, B = B-4r 3 p s p s , C = C + - 2 p s fi s , (15) 
a Z r° 



and A, B, and C are given by (|12p . Further, using the fact that the Weyl scalar has spin weight —2 [331 to expand 
^4 in terms of spin-weighted spherical harmonics 

*4 = X>, ro ftr)y_ 2 '' m (0,¥>). (16) 

in 1 

To deal with spin- weighted harmonics, it is useful to introduce the so called "eth" and "eth-bar" operators (42l. l43j] : 

( d d \ 
S s = - — + % csc 8 - s cot 8) = 9 + s cot 6», (17) 



88 dip 

( d d \ 

9 S = — I ttt — i csc 8 — — h s cot 8 = 9 — s cot 6, (18) 
\o8 dip J 

which have the property of raising and lowering the spin weight of the function on which they act. It is straightforward 
to show that the angular operator in (|11[) can be expressed in terms of the eth operators as 



(19) 
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The Y_2' m are the eigenfunctions of this angular operator: 

L 9v Y_ 2 l > m = B_i 9_ 2 F_ 2 <>" 1 = - (I - 1) (I + 2) y_ 2 <>"\ (20) 

In this way, the perturbation equation (|13[) transforms into a 1+1 equation for each (I, m)-mode of the gravitational 
perturbation, 

X-o-i) a + 2)] = 0, (21) 

where we have dropped the indices Z, m. In our numerical work, we have used a first order reduction of this equation 
as derived in the Appendix. 

III. INCLUDING NULL INFINITY IN THE NUMERICAL DOMAIN 

In the previous section we followed the standard procedure to calculate gravitational perturbations on a spherically 
symmetric background spacetime using the Teukolsky formalism with the additional feature that we allowed an 
arbitrary spherically symmetric foliation of the Schwarzschild background (see pij for the corresponding generalization 
in the Regge- Wheeler- Zerilli formalism). If we now insert the standard Schwarzschild gauge functions into equation 
pip we get the familiar Bardeen-Press equation. Our aim, however, is to include null infinity in the numerical domain. 
In this section we show how this can be done. 

A. Choice of a hyperboloidal foliation 

To include null infinity in the setting of a space-time splitting we first choose a hyperboloidal foliation by introducing 
a new time coordinate, r, whose level sets are hyperboloidal surfaces. Following [31j , we consider the class of spherically 
symmetric constant mean curvature (CMC) foliations of Schwarzschild spacetime (45l . |4(| . The transformation from 
the standard Schwarzschild time coordinate to the time coordinate of a spherically symmetric CMC-foliation can be 
written as r = t — h(r), where r is the Schwarzschild area radius, h{r) is the height function and t is the standard 
Schwarzschild time coordinate. The height function is not known in closed form. Its radial derivative is given by 



h'{r)= {l _T )p{r y - h ere J(r) = ^ - ±, and P(r) := ^ J{rf + (l - ^) . (22) 

Here, the mass of the Schwarzschild black hole is denoted by m. The foliation parameters are the mean extrinsic 
curvature K, and a constant of integration c. The global behavior of CMC-surfaces depends on the foliation parame- 
ters. We choose these parameters such that the surfaces of the foliation come from future null infinity, pass the event 
horizon above the bifurcation sphere and run into the future singularity as depicted in Figure[Tl[47|. The causal nature 
of the hyperboloidal foliation shown in the conformal diagram explains why we do not need a two stage calculation of 
the waveform in contrast to the characteristic case [13, EH ■ A further convenient feature of this foliation is that the 
natural coordinates are adapted to time symmetry in the sense that the timclikc Killing vector field of Schwarzschild 
spacetime dt is given by d T in the new coordinates. Therefore, timelike curves with constant spatial coordinates can 
be regarded as worldlines of natural observers at constant distances. 

In coordinates adapted to the CMC slicing, the standard Schwarzschild metric is obtained in the form 

/ 2m\ 2 2J(r) 1 ,2^2,2 

a = -( 1 -—) dT - p^y dTdr + j^) dr + r da > 

By comparing with (j4]) we see that the gauge functions read 

a(r)=P(r), (3(r) = —J(r) a(r), 7 (r) = — !— . (23) 

a[r) 

An important difference to Cauchy foliations should be noted at this point. The asymptotic behavior of an asymp- 
totically flat spacetime written in a Cauchy foliation implies a ~ 1 and (3 ~ as r — > oo. The same spacetime written 
in a hyperboloidal foliation implies a different asymptotic behavior for the gauge functions, namely 

a - 0(r) and [3 - 0(r 2 ) as r -> oo. (24) 



In the next subsection, we will sec that this behavior is important for the regularity of our equations. 



Figure 1: Penrose diagram of a CMC-foliation in 
Schwarzschild spacetime with m — 1/2, c = 1,K = 1. The 
dashed lines represent Killing observers. 



Figure 2: The function A^fO. 2 from (|28[) for the same hy- 
perboloidal foliation as in the diagram. We see that C/Q. 2 
vanishes at null infinity. 



B. Compactification 



To include the asymptotic ends of the surfaces of our hyperboloidal foliation, we introduce a compactifying coordi- 
nate p that is related to the standard Schwarzschild area-radius r via 

9 - 9 where f2 := 1 - p. (25) 



We choose the black-hole mass as m = 1/2, the horizon is thus located at p = 1/2. Future null infinity, is located 
at p = 1. Other choices are, of course, possible. Before writing the Bardeen-Press equation in the compactifying 
coordinate p, we factor out the singular asymptotic behavior of our hyperboloidal gauge functions (|2"4"|) in order to 
deal with regular functions. Considering the asymptotic behavior of the conformal factor, ~ 1/r, we define a := Q a 
and (3 := VL 2 j3. The rescaled gauge source functions become 



_2r^O\^ p = _ j5i where j :=nj= ^_£^!, 



These functions are obviously regular at {Q = 0}. Furthermore, they satisfy 



ft 



2 ^ )=0{n 2 ) with /3U+<0. (27) 



a 2 



We note that this behavior with the choice of the conformal factor (|25|) is specific to scri-fixing gauges and cannot be 
expected for all hyperboloidal foliations. It is related to the fact that the normal to level sets of r lies in ^ + [47] ]. 
The definitions of the auxiliary variables given in (| A. 1[) imply 

1 O 2 

= d r (j> = n 2 d p <j) = : fi 2 V>, and tt = -^{d t (f> - ftdM = -^(d t (j) - /3?/>) =: ^ 2 7f- 

or or 

We write the system (|A.2[) with respect to the functions <5, /3, ?/>, 7f by factoring an fi 2 term out of the equations. The 
resulting first order, symmetric hyperbolic system of equations reads 

d T <t> = a 2 7f + 

d T ij; = d p (a 2 7T + , 

d T TT = d p {a 2 Tp + /3tt) + A^tt + A^ip + - 

where A^,A^, and A^ are functions of p. It turns out that A n and A^ are conformally covariant. Writing A n and 
A^p from (|A.3[) in terms of the rescaled gauge functions and the compactifying coordinate p results in 



I 

p J ' \oi- 



Air = -4 (a 2 + P)dJ]n-\, A 4 , = A w - A a 2 3 P ^ 
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To show the regularity of our equations, we need to show that falls off with il 2 . We have by (|A.3[) 



+ ^ 02 „2 9 p\ -)■ ( 28 ) 



The explicit regularity of the above expression follows from the behavior of our rescaled gauge functions (|27|). specif- 
ically from a 2 + j3 ~ 0(f2 2 ) near null infinity. Indeed, as we can see in Figurc[2]the expression on the right hand side 
of (|28|) is not just regular but vanishes at null infinity for the foliation given in ([26]) . We note that scri- fixing plays an 
essential role here. 




IV. NUMERICAL STUDIES 



A. Initial data 



As we are not interested in the first (transient) part of the signal, the initial data can be chosen quite freely. We 
take a simple Gaussian as the initial perturbation and set 

0(O,p) = ae-^) 2 ^ 2 , ip(0,p) = d p <t>(0,p), tt(0,p) = 0. 

Here, p c is the center of the Gaussian pulse, a is its amplitude and a is its width. We set p c = 0.7, a = 1 and a = 0.05 
on the domain p G [0.495, 1]. The solid line in Figure [7] shows the initial data. We have tested that different choices 
of initial data do not have a relevant influence on the features we discuss in the following. Note that the vanishing of 
7r(0, r) does not mean that the data is time symmetric, because the initial surface is future hyperboloidal (spacelike 
while approaching future null infinity). We focus on the typically dominant I = 2 mode, corresponding to the lowest 
allowed angular momentum number for a spin-2 field because higher modes are weaker due to their faster decay rates 
and longer quasinormal mode ringing phases. 



B. Numerical algorithms 

We adopt the method of lines approach and discretize time and space separately, using 4th order Rungc-Kutta time 
integration and finite differencing with 4th, 6th and 8th order accurate stencils. We add Kreiss-Oliger type artificial 
dissipation to the evolution equation for n to suppress numerical high-frequency waves |48| . For a 2p — 2 accurate 
scheme we choose an operator (Q) of order 2p as 

h 2p-l 

where h is the grid size, D± are the forward and backward finite differencing operators and e is the dissipation 
parameter. Unless otherwise stated, we set e = 0.07. 

The inner boundary is a spacelike surface inside the event horizon where no boundary data is required. We excise 
the interior domain near the singularity. For the class of CMC surfaces we have chosen, excision is necessary as there 
is a minimal surface inside the event horizon where our coordinates break down. The outer boundary is at null infinity 
where one-sided finite differencing is applied. The one-sided finite differencing is of the same order as the interior one. 

The convergence of the code can be seen in Figure [3] For this plot, a three level convergence analysis has been 
performed with 505, 1010 and 2020 grid cells using finite difference operators of 4th, 6th and 8th order. For the 
6th and 8th order convergence tests, we used quadruple precision because the numerical error was below machine 
accuracy at late times when double precision was used. At very late times, convergence is lost due to accumulation 
of numerical errors (see the bending down of the curves in Figure [S] at late times) . 

Figure 3] shows the I/2-norm of the relative errors between the above mentioned simulations in different resolutions 
and finite differencing orders. Using an 8th order accurate finite differencing improves the final error by 8 orders 
of magnitude in the high resolution runs confirming that high accuracies can be achieved with high order finite 
differencing in combination with a relatively small amount of grid points. 

The simulation domain in the compactifying coordinate reads p £ [0.495,1]. The Courant factor Ai/Ap in the 
simulations has been chosen to be 4. We can choose a Courant factor that is larger than I because the coordinate 



8 




50 



100 



150 



200 



r/m 




Figure 3: Convergence in the Z/2-norm indicating 4th, 6th 
and 8th order convergence for the corresponding finite dif- 
ference operators. The convergence factor Q is calculated 

by Q = log 2 ^edZ^hiJ'n ■ 



Figure 4: i2-norm of the differences between numerical solu- 
tions obtained using different resolutions, i.e. \\cj> low — (f> rrled \\ 
and \\(f) med -(f) hlgh \\. Solid curves denote errors for 4th order 
finite difference operators, pointed curves for 6th order and 
dashed curves for 8th order. 




Figure 5: The dashed and the solid curves denote coordinate speeds of in- and outgoing characteristics respectively, for the 
foliation characterized by (|26p . The horizon is located at p = 0.5. We see that the characteristic speeds are less than unity 
which allows us to choose a large Courant factor. We observe that there is little variation of the speeds over the grid which has 
infinite physical extent. 



speeds of characteristics are less than 1. For outgoing and ingoing characteristics, the characteristic speeds with 
respect to our foliation are v± = ±a 2 ~ (3. Figure [5] shows the characteristic speeds along the grid for K = 1. The 
allowed Courant factor depends on the value of K as K is related to the coordinate speed of characteristics. For 
example, the coordinate speed of outgoing characteristics at reads in our foliation = 2K 2 /9. Therefore, 

a large value for K implies a small Courant factor [47| ■ We plot in Figure \E\ the speeds for K = 1 and see that 
v± < 1 which allows At/Ap > 1. We also see that the ingoing characteristic speed vanishes at as expected, 
correspondingly we do not need to impose outer boundary conditions. Note also that there is relatively little variation 
of the characteristic speeds over this grid which has infinite physical extent. 
In the studies presented below, we choose m = 1/2, c = 1 and K = 1. 



C. Quasinormal mode ringing 

Figure [5] shows the gravitational perturbations of a Schwarzschild black hole in a double logarithmic plot with 
respect to three observers located at r « {2.5m, 18m} and at J? + . We see that, as expected, the quasinormal mode 
frequencies arc the same for each observer. It docs not matter for the calculation of these frequencies at what distance 
from the black hole the perturbations are measured. The end of the quasinormal ringing phase, however, is observer- 
dependent. The polynomially decaying part of the solution starts earlier for an observer at ^ + than for a nearby 
observer. This suggests that for the calculation of quasi-normal mode frequencies, one may get better results by using 
the waveform as measured by a nearby observer. 

We expect the quasinormal mode ringing to have the form 

4> p {t) = a p e- U2T sm(u 1 T + <p p ). (29) 

Here, uj\ and 0J2 are the mode frequencies, a p is the amplitude and <j> p is the phase of the wave signal. In contrast to 
the mode frequencies, the amplitude and the phase do depend in general on the observer's location. 
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Figure 6: Quasi normal mode ringing of Schwarzschild 
spacetime excited by a Gaussian gravitational perturbation 
with 1 = 2 measured by three observers. Looking to the 
tail part of the signal, the observers are located from top to 
bottom at ,y + , r ~ 18m and r w 2.5m. We see that the 
ringing phase ends at different times for different observers 
and goes into a polynomial decay. 



4> 




Figure 7: The numerical solution on the grid during 
the quasinormal ringing phase plotted at different times 
and rescaled for visibility. The solid curve is just the 
initial data. The dashed curves are from top to bottom: 
{10 x </}(80m, p), 100 x 0(lO6m, p), 1000 x 0(122m, p)}. 
We see that the wave package is distributed quite uniformly 
across the grid. 



To demonstrate the efficiency of our method, we calculate the quasinormal mode frequencies during the ringing 
phase using various resolutions. The frequencies are compared to the result from Leaver's continued fraction method 
[49j ]. For our choice of black hole mass of m — 1/2 they read (wi,^) = (0.747343,0.177925). Tabic U shows the 
measured frequencies for various grid resolutions for 4th and 6th order finite differencing operators along with the 
relative errors. The relative error for w° um ' 2p ' in a numerical simulation with spatial finite differencing of order 

2p is calculated by Su)^ p ^ = Icj* 1 " 111 ^ 2 ^ — u>i\/u>i. The fitting is performed using a simple least squares method on the 
interval t <G [80m, 180m] with an initial guess of the frequencies. 



Cells 




w num(4) 


<M 4) 


s4 4) 


w num(6) 


w num(6) 


5<j[*> 


s4 6) 


25 a 


0.825455 


0.157916 


10~ 2 


1.1 x 10~ 2 


0.777778 


0.168439 


4.1 x 10~ 2 


5.3 x 10~ 2 


50 


0.746988 


0.177897 


4.7 x 10~ 4 


1.5 x 10" 4 


0.747341 


0.177917 


2.3 x 10~ 6 


4.4 x 10" 5 


100 


0.747325 


0.177927 


2.4 x 10~ 5 


1.1 x 10~ 5 


0.747344 


0.177923 


1.9 x 10~ 6 


9.5 x 10~ 6 


200 


0.747343 


0.177924 


2.6 x 10~ 7 


4.9 x 10" 6 


0.747344 


0.177924 


1.8 x 10~ 6 


5.4 x 10" 6 



"For the runs with 25 grid cells, a relatively high dissipation parameter of e = 0.2 has been used. 

TABLE I: Numerical quasinormal mode frequencies Ii ;" um ' 2p ' i n various resolutions measured at r ~ 2.5m and the relative 
errors <5w^ 2p ' 1 calculated in a numerical simulation with spatial finite difference operators of the order 2p. We perform a simple 
least square fit to (|29p with initial guess. We see that even with low resolutions, we are able to measure the quasinormal-mode 
frequencies with a high accuracy indicating the efficiency of our method. 

We see that already with 50 cells we are able to calculate the quasinormal mode frequencies with a relative error 
of about 10~ 4 using 4th order finite differences and 10~ 6 using 6th order finite differences. The use of higher order 
accurate finite differencing does not increase the computational cost significantly. We observed only about a 5% 
increase in the run time for the 6th order finite differencing compared to the 4th order one. We note that for a 
relative error at the order of 10 -6 , the accuracy in the frequency is rather dominated by the fit than by numerics. 
This explains, in part, why the relative error in the simulation with 200 grid cells seems smaller in the 4th order case 
than in the 6th order case. Using 8th order finite differencing gives similar results as the 6th order one and does not 
increase the accuracy in the frequencies considerably. 

The fact that using a very small number of cells is already sufficient for such accuracy in the frequencies demonstrates 
the efficiency of the hyperboloidal approach. An explanation for this efficiency can be found in the way the solution 
looks like. Figure [7] shows the solution at three different times during the quasinormal ringing phase. We see that 
the initially localized perturbation is distributed uniformly on the grid. This uniform distribution implies that there 
is no need to use a high number of grid points to resolve the outgoing wave package. For accuracy, it is more efficient 
to use a high order scheme with a small number of grid points in this case. In evolutions with Cauchy-type foliations, 
however, the outgoing wave package is a long wave-train with the typical quasinormal ringing form that requires good 
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resolution. We note that the observed broadening in the hypcrboloidal approach is stronger for larger values of K. 
For smaller values we get a similar localization like when Cauchy-type foliations are used. This is because K can be 
regarded as a measure of how close the behavior of hyperboloidal surfaces is to being characteristic (42| ■ 



D. Polynomial decay 



Figure [8] shows the polynomially decaying part of the solution starting from r = 500m as seen at future null infinity, 
at r « 60m and at the future horizon, TC + , in an evolution with 4040 grid cells. We observe that the scattered 
perturbations are stronger for the observer at and also the slope of the curves are different, indicating that the 
decay rate for a far away observer is less than the decay rate for a close by observer. An intuitive explanation for this 
difference is given by the picture that the tail is generated by backscattering off the nonvanishing curvature of the 
background spacetime. An observer who is farther away sees more scattering than a nearby observer. Therefore, the 
amount of scattered waves measured by a nearby observer is less than for a farther one, consequently the scattered 
part of the signal is "stronger" for the far away observer (modulo the natural fall-off) and decays slower. The natural 
fall-off is not seen in Figure [5] as the variable we plot corresponds to the real part of 'J 4 rescaled with r. 




Figure 8: The polynomially decaying part of the solution 
in a double-logarithmic plot starting from r = 500m. The 
upper curve gives the solution at null infinity, the middle 
one at r ~ 60m and the lower one at the future horizon 7i + . 
It can be seen that the decay of the solution along J ;+ is 
slower than near the black hole. 



Figure 9: The local power index p p (t) for observers located 
from top to bottom at {J? + , 800m, 260m, 110m, 50m, 25m}. 
We see that after about r = 4500m the local power index for 
nearby observers bends down. This behavior is seen earlier 
for lower resolutions and is only a numerical effect indicating 
loss of convergence due to loss of accuracy. 



We see from Figures [5] and [5] that we are able to follow the solution for more than 20 orders of magnitude and 
calculate the predicted decay rate with high accuracy. To achieve such accuracy, we use 8th order finite difference 
operators. The exponentially decaying (quasinormal mode ringing) part of the solution is calculated in quadruple 
precision. The reason for choosing 8th order accurate finite differencing is due to the very small effects that we are 
trying to measure. Note that for any given finite difference order, there are ghost potentials which may be mistaken 
as tail solutions We have tested that the tail seen in Figure [5] is not a numerical artefact by comparing solutions 
in various resolutions. 

An accurate picture of the decay rates can be obtained from Figure [9l Here we plot, for various observers, the local 
power index defined as 

amr 

The function p p (r) becomes asymptotically the exponent of the polynomial decay of the solution. At late times, 
our numerically measured decay rate becomes —6.01 for an observer at ^ + and —7.002 for an observer near i + , in 
accordance with the well-known prediction of 2~ (2l+2 ^ at y + and of 2~ (2l+3 ^ near i+ |5l| (see also [H, [H, [H, H3, [H, 

MM). 
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V. DISCUSSION 



In this paper we have revisited the problem of scattering pulses of gravitational waves off a Schwarzschild black hole 
with the aim of getting guidance in the development of a hypcrboloidal approach to numerical relativity. We have 
used constant mean curvature slices on the background of a Schwarzschild spacetime (4f| EH in a scri- fixing gauge [47} 
to numerically solve a hyperboloidal initial value problem (3lj for the Bardeen-Press equation. The numerical study 
of gravitational perturbations at null infinity follows, in some sense, the spirit of [13, [HHHI w ho studied compactified 
null foliations. 

Motivation for taking a hyperboloidal approach comes from the well-known shortcomings of existing methods based 
on Cauchy-type and characteristic foliations for this problem. The drawbacks of using Cauchy-type foliations can be 
summarized as follows: (i) waste of computational resources when one is interested in the waveform as seen by far 
away observers, (ii) spurious reflections from the outer boundary contaminating the waveform in long time evolutions 
due to unknown boundary data, (iii) not being able to read off the waveform at null infinity. Using a characteristic 
foliation solves these problems, but leads to new difficulties: (i) the causal relationship between the surfaces seems to 
require a two-stage approach [Tt], EH , (ii) choice of variables for the numerical simulations is not straightforward due 
to loss of resolution, and most importantly, (iii) extension beyond spherical symmetry is very complicated. Especially 
in highly dynamical situations where formation of caustics is expected [58j , it is not yet clear how the characteristic 
approach can be implemented [59| . 

The hyperboloidal approach hopes to combine the best properties of currently used approaches in the field, namely 
the flexibility of Cauchy-type surfaces and the asymptotic behavior of characteristic surfaces, and avoids their problems 
- however it is also much less developed so far. 

We believe that the problem of treating gravitational perturbations on a fixed spherically symmetric background, 
while rather simple, is an important step in the development of algorithms for the hyperboloidal initial value problem. 
We have chosen to solve the Bardeen-Press equation describing perturbations of Schwarzschild spacetime in terms of 
the Wcyl component ip4, partly because this naturally generalizes to the much more interesting axisymmetric case via 
the Teukolsky equation. 

We demonstrated in this paper that compactification along hyperboloidal surfaces using the Teukolsky formalism 
on the background of a Schwarzschild spacetime leads to regular equations. This regular behavior can be expected 
on the background of any weakly asymptotically simple spacetime (in the sense of Penrose [60l l6l|). The resulting 
equations can be solved numerically to study quasinormal mode ringing very efficiently. The tail part can be calculated 
accurately with sufficient numerical resolution. There are no reflections from the numerical outer boundary. The only 
approximation we use to calculate the gravitational waveform beside the linearization is the numerical discretization 
of the equations. The observed uniform distribution of outgoing wave packages along the grid and the numerical 
experiments we performed suggest that the use of high order methods to solve the linear perturbation equations lead 
to an efficient code that requires a moderate number of grid points while achieving good accuracy in the results. 

We want to conclude with a brief discussion of open problems for future work. It may be interesting to compare 
the approach for curvature-based and metric perturbations using the Chandrasckhar transformation along the lines 



of 14{ |. Current work is under way to study metric perturbations using the Regge-Wheeler-Zerilli formalism in a 



hyperboloidal setting [62|. The hyperboloidal approach to study gravitational perturbations at null infinity should 
generalize rather straightforwardly to Kerr spacetime, for example by using the foliations presented in |47| . 

The difficulty of generalization to the fully nonlinear Einstein equations remains unclear. In contrast to the test 
field case, the conformally rescalcd Einstein equations are not explicitly regular. There, a scri-fixing gauge can be 
constructed in a hypcrboloidal initial value problem for the generalized harmonic reduction of the Einstein equations 
such that the formally singular terms attain regular limits at null infinity [63j . The open question is whether such 
limits can be calculated numerically in a stable way. The asymptotic behavior of the gauge functions ([27)) need to be 
taken into account if a regularization along the lines of [63[ needs to be done for the ADM or BSSN systems based 
on scri-fixing. A very promising alternative is to adopt an elliptic-hyperbolic formulation of the problem, and apply 
boundary conditions for gauge quantities as an elliptic Dirichlet problem, as suggested in [37] ]. 
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APPENDIX: SYMMETRIC HYPERBOLIC FORM OF THE BARDEEN-PRESS EQUATION 

We solve the Bardeen-Press equation in first order symmetric hyperbolic form using a spherically symmetric CMC- 
foliation of Schwarzschild spacetime with the gauge functions given in (|2"5|) . The relation a = I/7 simplifies the 
calculation, so we will assume it in the following. In such a gauge, the spinor coefficients take the form 

a 2 - f3 a 2 +f3 
2 a z r r 

To bring the radial part of the Bardeen-Press equation (|21[) to first order form, we introduce the auxiliary variables 

if) := d r <f>, and n := -^-(<9 t — (3d r <ft). (A.l) 
or 



Then the following linear, symmetric hyperbolic system of evolution equations is obtained from (|14[) . (|15p 

d t <t> = a 2 7t + /3ip, 

d t vb = d r (a 2 Tr + ^), (A.2) 
d t n = d r (a 2 ip + f3n) + A„ tt + A, 4 , ip + (a^ - Aj fa 

with A = (I - 1) (I + 2) and 



A„ = -4 (c^ + /3)cU In - 



r 



A^ = A^-Aa'dr I 4 ) . ( A -3) 



r \ or / r \ \ or / \ a J \ r / J \ a + p 
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